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The morphology of presynaptic specializations can vary greatly ranging from classical 
single-release-site boutons in the central nervous system to boutons of various sizes 
harboring multiple vesicle release sites. Multi-release-site boutons can be found in 
several neural contexts, for example at the neuromuscular junction (NMJ) of body wall 
muscles of Drosophila larvae. These NMJs are built by two motor neurons forming two 
types of glutamatergic multi-release-site boutons with two typical diameters. However, 
it is unknown why these distinct nerve terminal configurations are used on the same 
postsynaptic muscle fiber. To systematically dissect the biophysical properties of these 
boutons we developed a full three-dimensional model of such boutons, their release 
sites and transmitter-harboring vesicles and analyzed the local vesicle dynamics of various 
configurations during stimulation. Here we show that the rate of transmission of a bouton 
is primarily limited by diffusion-based vesicle movements and that the probability of vesicle 
release and the size of a bouton affect bouton-performance in distinct temporal domains 
allowing for an optimal transmission of the neural signals at different time scales. A 
comparison of our in silico simulations with in vivo recordings of the natural motor pattern 
of both neurons revealed that the bouton properties resemble a well-tuned cooperation 
of the parameters release probability and bouton size, enabling a reliable transmission 
of the prevailing firing-pattern at diffusion-limited boutons. Our findings indicate that the 
prevailing firing-pattern of a neuron may determine the physiological and morphological 
parameters required for its synaptic terminals. 



Keywords: neuromuscular junction, boutons, modeling and simulation, structure-function relationships, 
morphology, firing pattern, adaption 



1. INTRODUCTION 

Chemical synapses are specialized compartments of neurons 
for neuron-to-neuron or neuron-to-muscle communication that 
exist in various morphological layouts (Rollenhagen and Lubke, 
2006) ranging from classical single-release site synapses at cortical 
neurons (Gulyas, 1993) over intermediate-sized multi-release- 
site boutons of the calyx of Heldt (Satzler et al., 2002) or 
the arthropod neuromuscular junctions (NMJs) (Bradacs et al., 
1997) to giant multi-release-site terminals at vertebrate NMJs 
(Dreyer et al., 1973). While the principal processes underly- 
ing synaptic transmission and plasticity at these synapses are 
rather well understood, the role that the various morpholog- 
ical layouts may play in shaping synaptic functions is cur- 
rently unclear. One reason for this astonishing discrepancy 
may be the difficulties to systematically control the morpho- 
logical and physiological configuration of a given synaptic 
system. 



Prominent examples for two highly related but morphologi- 
cally distinct multi-release-site boutons are the well-characterized 
glutamatergic synapses at larval NMJs of Drosophila (Jan and Jan, 
1976; Atwood et al, 1993; Schuster, 2006; Collins and DiAntonio, 
2007). Larval body wall muscles of Drosophila are typically inner- 
vated by two motor neurons of which one forms large spherical 
and somewhat variable type lb boutons and a second motor neu- 
ron forms smaller and more regular type Is boutons. Both types 
of boutons harbor multiple glutamatergic release sites that are 
spaced according to a strict nearest-neighborhood relationship 
(Atwood et al, 1993; Meinertzhagen et al, 1998; Sigrist et al, 
2003) and that show discrete differences in their release prob- 
abilities (Kurdyak et al., 1994; Cooper et al, 1995). In spite of 
their close functional relationship, it has so far not been possible 
to systematically assess why both motor neurons maintain their 
obvious morphological and physiological differences and what 
their relevance is for the animal's behavior. 
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FIGURE 1 | Functional reconstruction of synaptic boutons of 
Drosophila in silico. (A) Confocal image of a Drosophila larval 
neuromuscular junction (muscle 13) showing the synaptic terminals of two 
motor neurons that form type Is or lb boutons, respectively. Release sites 
(green) were labeled with an antibody recognizing the active zone protein 
Bruchpilot (Wagh et al., 2006). The axonal surface (red) was labeled by 
anti-HRP immunoreactivity. Right panel: magnification of representative Is 
and lb boutons (arrow and arrowhead) displaying an average diameter of 2 
and 3\im, respectively. Scale bars: 5|im (left) 1.5 jun (right). (B-D) Release 
sites (green) and the vesicular diffusion area (red) are represented in silico 
by tetrahedral volume grids (Bastian et al., 2000). Release sites are 
positioned according to nearest neighborhood laws derived from 
ultrastructural data (Meinertzhagen et al., 1998; Sigrist et al., 2003). The 
bouton center is free of vesicles (Denker et al., 2009). (E,F) Simulation of 
vesicle dynamics in a lb bouton (3 \i.m) during 30 Hz or 60 Hz stimulation 
(P 0 = 7%) with initial vesicle density of 275 |im~ 3 . Encoded in color is the 
spatially resolved vesicle density (in vesicles/iim 3 ) at the indicated stimulus 
number (599-601 ). Note that 60 Hz-stimulation rapidly depletes most 
release sites and their immediate surroundings of releasable vesicles 
(white) while more distant areas maintain high vesicle densities (red). 



In order to shed light onto the poorly understood structure- 
function relationship of such discrete bouton-configurations, we 
used a combined theoretical and experimental approach to dissect 
the parameters that are biophysically important to support the 
natural function of discrete bouton configurations. We therefore 
developed and verified a three-dimensional model of synaptic 
bouton functions of larval NMJs of Drosophila and used it to 
systematically assess the meaning of experimentally rather inac- 
cessible bouton parameters like bouton size, vesicle diffusion and 
vesicle release probability. We found that these parameters affect 
bouton output and its reliability of transmission on different time 
scales and that their individual combination in real boutons is 
predictive of the prevailing firing pattern that a given bouton type 
experiences. 

2. RESULTS 

2.1. MORPHOLOGICAL AND PHYSIOLOGICAL PARAMETERS TO 
GENERATE A MODEL OF GLUTAMATERGIC BOUTONS 

To systematically assess the influence of morphological and phys- 
iological constraints on the performance of multi-release-site 
boutons we first established a functional reconstruction of such 
boutons in silico that aimed at simulating the vesicle dynam- 
ics within a bouton during trains of action potentials. From 
confocal image stacks of NMJs of third instar Drosophila lar- 
vae (Figure 1A) and representative measurements of bouton 
dimensions of type Is and type lb synaptic terminals, we first 
estimated the typical bouton diameter of type Is and type lb 
boutons to 2 and 3 u.m, respectively (right panel in Figure 1A). 
In addition, some terminal boutons showed diameters of up to 
5 |xm. These different bouton morphologies were constructed 
in ProMesh (Reiter and Wittum, 2014), and discretized with 
a tetrahedral volume grid. The geometries represented spheres 
with a volume spared out at the center of each sphere repre- 
senting organelle occupation, hence an area not occupied by 
vesicles (Figures 1C,D, 2). Each bouton geometry consisted of 
a zone Q,q where vesicles were motile as well as exocytosis 
zones £2i . . . £2 n where vesicles are motile and could be released 
(Figures 1B-D, 2). Taking into account the total number of vesi- 
cles that can be released from a NMJ during vesicle depletion 
experiments (Delgado et al, 2000) and data from ultrastructural 
analyses (Sigrist et al, 2002) we estimated a concentration of 
150-500 vesicles per |xm 3 within a bouton at rest (see Materials 
and Methods for details). The positioning of the vesicle release 
sites on the surface of these spherical boutons was guided by serial 
ultrastructural reconstructions of larval boutons from which we 
determined the densities of T-bar harboring active zones between 
0.37 and 0.46 |xm 2 (Meinertzhagen et al, 1998; Sigrist et al, 
2002, 2003) with an average active zone diameter of 300-400 nm 
(Meinertzhagen et al., 1998; Sigrist et al., 2002, 2003). Based on 
these data we assumed 7 equally spaced active zones per Is bou- 
ton, 10 release sites per lb bouton and up to 14 release sites in a 
large terminal bouton. 

Synaptic active zones were represented by a multiedge- 
surface approximating a cylindrical structure with a diameter of 
300-400 nm and a volumetric depth of 200 nm, which defines 
a volume that harbored roughly 6-8 vesicles, depending on the 
starting vesicle densities, which were chosen between 150 and 



500 vesicles per u,m 3 . An overview of the parameter ranges used 
in the simulations is given in Table 1. Where simulation data 
and results are shown, we specify which values within the given 
ranges were chosen. Within the active zones we made a tran- 
sition from vesicle motion at the continuous scale to discrete 
vesicles that undergo exocytosis (see Materials and Methods for a 
detailed description). At this stage of the model we did not incor- 
porate voltage-gated calcium channels on the plasma membrane 
and presynaptic calcium dynamics and its role in triggering vesi- 
cle fusion but instead modeled action potential triggered release 
of vesicles by the term / syn that described a stochastic vesicle 
release process controlled by the stimulation frequency on and 
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FIGURE 2 | Two-dimensional sketch of the bouton model. Region S2o 
denotes the bouton area in which vesicle motion is governed by diffusion. A 
cutout region in the center represents space occupied by organelles, such 
as mitochondria. Q-i . . . Q 4 define synaptic active zones. In these 
subdomains, vesicles diffuse and are subject to exocytosis modeled by the 
term f SY „, that also includes the transition from the continuous to discrete 
scale in £2i . . . £2,}. The organelle as well as bouton boundaries are treated 
as no-flux boundary conditions. 



Table 1 | Overview over the bouton configurations and parameter 
values used in this study. 



Boutons 


Is 


lb 


Terminal bouton 


Bouton diameter 


2 u,m 


3 u.m 


5 u,m 


Number of synapses 


7 


10 


14 


Cutout Radius 


0.5 u.m 


0.8 [im 


1.3 urn 


Values 


Range 


Stepwidth 




Po 


1-90% 


1% 




ft) 


10-100 Hz 


5 Hz 




Vesicle density 


150-500 iim" 3 







a constant output probability P 0 throughout a given simula- 
tion. For literature-based simulations, we used P 0 = 29 and 7 
% assigned to Is and lb boutons respectively, values taken from 
Kurdyak et al. (1994); Cooper et al. (1995) instead of regulated 
P 0 values that are known to be involved in synaptic facilitation 
and depression. This restriction enabled us to dissect the role of 
each of these parameters on presynaptic vesicle dynamics under 
otherwise constant conditions. Although vesicle recycling has an 
effect on the transmission efficacy, we did not incorporate vesicle 
replenishment into the model in order to avoid introducing still 
unknown variables. Since the experiments also focus exclusively 
on preexisting vesicles, model and experiments are comparable. 



To implement into our model a meaningful simulation of 
vesicular motilities during trains of high-frequency action poten- 
tial firing we made the following considerations: at Drosophila 
NMJs presynaptic vesicles have been assigned to three function- 
ally distinct vesicle pools, the readily releasable pool, the recycling 
pool and the reserve pool of vesicles (Kuromi and Kidokoro, 2000; 
Denker et al., 2009). Although the number of vesicles belong- 
ing to each vesicle pool and their respective motilities are very 
different at rest (Kuromi and Kidokoro, 2000), stimulation fre- 
quencies above 10 Hz have been shown to be sufficient to recruit 
vesicles from all pools for release (Kuromi and Kidokoro, 2000; 
Wagh et al., 2006; Denker et al, 2009). Since the natural fir- 
ing pattern of larval motor neurons uses typically frequencies 
higher than 10 Hz (Figure 6) all vesicles within a bouton have 
to be considered mobile during natural activity patterns in vivo. 
It is, however, currently not clear whether vesicle motility dur- 
ing high frequency stimulation follows a simple diffusion process 
(Gaffield et al., 2006) or whether it underlies active transport 
mechanisms requiring cytoskeletal elements (Ryan, 1999; Jordan 
et al, 2005; Tokuoka and Goda, 2006; Seabrooke et al, 2010; 
Kisiel et al., 201 1). We have previously shown that the presynap- 
tic actin-based cytoskeleton within Drosophila boutons is highly 
dynamic with major rearrangements occurring on a rapid time 
scale (Steinert et al., 2006). These cytoskeletal rearrangements 
may result in a constant intermixing of presynaptic vesicles and 
hence in a dynamic and overall rather homogeneous distribu- 
tion of vesicles within the presynaptic terminal. Based on these 
and the above observations we first assumed that all presynap- 
tic vesicles are mobile during stimulation frequencies higher than 
10 Hz (Kuromi and Kidokoro, 2000). We, secondly, assumed that 
during such stimulation periods vesicle motility may be approxi- 
mated by undirected diffusion-like movements (Shtrahman et al., 
2005; Darcy et al., 2006; Fernandez-Alfonso and Ryan, 2008; 
Westphal et al., 2008). We therefore used an effective diffusion 
coefficient of D = 5 • 10~ 3 (im 2 /s (Shtrahman et al., 2005) that 
included stick-and-diffuse processes that have been shown to 
occur during interactions of vesicles with the cytoskelleton (Darcy 
et al, 2006; Fernandez-Alfonso and Ryan, 2008; Westphal et al., 
2008). 

The above considerations formed the basis for the formulation 
of the diffusion equation and the sink term that were solved on 
complex morphologies in three-dimensional space using numeri- 
cal discretization methods and solvers (see Materials and Methods 
section for details). All numerical methods were part of the simu- 
lation platform uG (Bastian et al, 2000). Numerical simulations 
of our model generated three-dimensional density profiles of 
presynaptic vesicles that allowed us to monitor the local distri- 
bution and dynamics of mature vesicles as a function of bouton 
diameter, the applied stimulation pattern and the vesicle release 
probability (Figures 1E,F, Supplemental Video SI). 

2.2. DISPR0P0RTI0NAL DEPLETION OF VESICLE RELEASE AT HIGH 
STIMULATION FREQUENCIES IS DUE TO LIMITED VESICLE 
DIFFUSION 

To verify the validity of the parameters of our bouton model 
we wished to compare its behavior in a task that is similar to a 
simple in vivo experiment. We therefore chose an experimental 
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FIGURE 3 | Disproportional depletion of eEJPs is due to limited 
availability of presynaptic vesicles. (A) Experimental design: larval filet 
preparations were pre-incubated in normal saline (NS) or NS containing 
2 \M Bafilomycin A1 (Baf A1 ), a blocker of the vesicular proton pump, to 
inhibit transmitter filling in recycling vesicles. Motor nerves were then 
continuously stimulated with a suction electrode at various frequencies 
(HFS) while we simultaneously recorded the membrane potential of the 
postsynaptic muscle fiber 6. (B,C) Representative traces of these 
recordings in the absence or presence of BafA1 taken at the indicated 
time-points of stimulation (stars: stimulation artifacts). (D) Quantification of 
eEJP-depletion dynamics during continuous 30 Hz stimulation with and 
without BafA1 incubation. Note, that eEJPs could be detected over at least 
15 min of 30 Hz-stimulation in controls without BafA1-treatment (B and 
circles) whereas eEJPs of BafA1-treated preparations depleted within 2 min 
(C and squares). (E-G) eEJP-depletion dynamics in BafA1-treated 
preparations caused by the indicated stimulation frequencies and 
expressed as a function of the stimulus number. (E) Shows the 
eEJP-depletion dynamics during continuous long-term stimulation. In (F,G) 
stimulation trains of 100 stimuli were separated by pausing intervals of 5s. 
Note that the depletion of eEJP-amplitudes was disproportionally faster 
during 60 and 80 Hz stimulations compared to 15 and 30 Hz stimulations. 
Short pausing intervals resulted in an almost complete recovery of 
eEJP-amplitudes at the beginning of the experiment (F). In later phases of 
the experiment (G) the eEJP recovery was significantly higher in 80 Hz 
stimulated NMJs than in 15 Hz stimulated NMJs (p<0.033 at stimulus 1800 
and p<0.019 at stimulus 1900). Data in (D-G) represent means ± SEM of 4 
to 6 independent preparations. 



vesicle depletion paradigm that is triggered by continuous high- 
frequency stimulation of larval motor nerves in the presence of 
a blocker of the vesicular proton pump, Bafilomycin Al (Kuromi 
and Kidokoro, 2000) (BafAl in Figure 3A). Synapses in BafAl- 
treated larval filet preparations fail to fill recycling presynaptic 
vesicles with the neurotransmitter glutamate and hence can be 
used in intracellular recording experiments of postsynaptic mem- 
brane potentials to exclusively visualize the release of preexisting 
mature vesicles. As expected, the amplitudes of evoked excitatory 
junctional potentials (eEJPs) of BafAl -treated larval preparations 
depleted significantly faster during continuous 30 Hz nerve stim- 
ulation than in untreated preparations (Figures 3B-D). The here 
observed rate of eEJP-deletion in BafAl -treated preparations was 
similar to the rate of vesicle depletion in temperature-sensitive 
mutants of dynamin (shibire ts ) (Siddiqi and Benzer, 1976) that at 
restrictive temperatures inhibit all endocytic processes (Delgado 
et al., 2000) indicating that in spite of the different molecular 
mechanisms of both treatments the measurable outcome was 
comparable. In addition, we observed that the eEJP-depletion 
progressed at rather slow rates during stimulation frequen- 
cies between 15 and 30 Hz, whereas at higher stimulation fre- 
quencies (60-80 Hz) the initial depletion was disproportionately 
faster (arrows in Figure 3E). This unexpected phenomenon was 
unlikely due to desensitization of postsynaptic glutamate recep- 
tors, since the quantal amplitudes (miniature excitatory junc- 
tional potentials, mEJPs) were unaltered under all tested condi- 
tions (data not shown). Instead, the latter observation suggested 
that a presynaptic limitation was responsible for the dispro- 
portionally rapid decay of eEJP amplitudes at high stimulation 
frequencies. 

An analysis of NMJs that transmitted naturally produced 
motor patterns in the absence of BafAl (Supplemental Figures 
S1A,D) revealed that these synapses were capable of maintaining 
their eEJP amplitudes throughout entire patterns (Supplemental 
Figures S1C,F) even when the firing frequencies were in the 
range of 60-80 Hz (Supplemental Figures S1B,E). This result 
demonstrated that high-frequency firing patterns could be faith- 
fully transmitted by the glutamatergic boutons of larval NMJs. 
It further suggested that the disproportionally rapid decay of 
eEJP amplitudes seen in BafAl -treated synapses might be due 
to limited availability of preexisting mature vesicles at active 
zones. If such a limited availability of transmitter filled vesi- 
cles would cause the disproportional synaptic depression during 
60-80 Hz stimulation rates one would expect that short paus- 
ing intervals between stimulation trains allow distant mature 
vesicles to functionally reoccupy release sites over time result- 
ing in recovered eEJP amplitudes. If BafAl -treatment would 
potentially unspecifically block other synaptically relevant pro- 
cesses such as calcium export or cytoskeletal organization that 
would disproportionally stronger affect high frequency trans- 
mission, eEJP amplitudes should continue to show depression 
after a paused stimulation. This assumption is contradicted by 
our results in Figures 3F,G showing transient depression at all 
stimulation frequencies excluding a persistent pleiotropic defect 
induced by the procedure. In experiments in which trains of 
100 stimuli were separated by pausing intervals of 5 s we found 
an almost complete recovery of the eEJP amplitudes at all used 



stimulation frequencies (15-80 Hz) immediately after the pause 
(Figure 3F). Furthermore, the longer this mode of stimulation 
continued the more the recovered eEJP amplitudes deviated 
between high and low frequencies (Figure 3G): eEJP recovery 
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was significantly larger in boutons after 1800 stimuli at 80 Hz 
than at 15 Hz (p < 0.033, n = 4) or after 1900 stimuli (p < 0.019, 
n = 4). Both observations were consistent with the idea that vesi- 
cle motility was limiting the availability of mature vesicles at 
active zones. This leads during high frequency stimulation to a 
rapid failure to reoccupy empty active zones with mature vesi- 
cles and thus to a rather small total number of released vesicles 
per train. The higher number of remaining mature vesicles is 
then responsible for the strong eEJP recovery even after pro- 
longed stimulation periods (triangle in Figure 3G). Conversely, 
during lower stimulation frequencies vesicle release and active 
zone reoccupation appear to be more balanced resulting in 
a higher total vesicle output per train and consequently in 



a weaker eEJP recovery after prolonged stimulation (circle in 
Figure 3G). 

2.2. 1. Verification of the model 

To further verify the quality of our model we systematically com- 
pared simulated and experimental vesicle depletion curves of 
different conditions (Figure 4). For the model we used as before 
parameter values from previous experimental studies (Kurdyak 
et al, 1994) and set up our model with 6 lb and 13 Is boutons 
with P 0 values of 7 and 29% for the release sites of type lb and 
type Is boutons, respectively (Kurdyak et al., 1994). Under these 
conditions the in silico simulations of the vesicle depletion exper- 
iments approximated the experimental observations in that the 
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FIGURE 4 | Model validation at different frequencies and stimulation 
protocols. (A) Comparison between experimental and simulated data at the 
frequencies 15, 30, 60, and 80 Hz. In order to compare simulation and 
experiment, the initial vesicle density was computed from the experimentally 
recorded vesicle release for each frequency and thus initial vesicle densities 
were set to 500, 280, 190, and 150 vesicles/(i.m 3 respectively in the 
simulations. Note that in the beginning 1000 stimuli at 15 Hz the depletion is 
considerably faster in real boutons than in the simulation. This may indicate 
that 15 Hz stimulation is initially not sufficient to efficiently mobilize all 
presynaptic vesicles e.g., from the reserve pool of vesicles. This limitation is 



overcome at higher stimulation frequencies. (B-E) Experimental (dotted) and 
simulation (lines) data in stop and go stimulation protocol, with 100 stimuli 
and 5 s pause, show decline in in the postsynaptic amplitudes (B,C: 
simulation data averaged over 25 stimuli; D,E: simulation data averaged over 
5 stimuli). The plots show that with increasing frequency (from B-E) 
amplitudes decrease slower, due to higher exocytosis failure rates at high 
frequencies caused by diffusion limited active zone replenishment. Note that 
similar to (A) the experimental depletion is faster than the in silico depletion. 
As above this may be due to a not efficient mobilization of presynaptic 
vesicles at 15 Hz. 
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simulated total amount of vesicular release decreased at various 
stimulation frequencies with similar time courses as in the in vivo 
experiments (Figure 4). Intriguingly, we identified the same two 
features that together could explain the disproportionally faster 
depletion of vesicle release at high stimulation frequencies with 
our model simulation of the local vesicle dynamics within a bou- 
ton (Figures 1E,F): first, 600 stimuli applied at 60 Hz depleted 
most release sites and their immediate surroundings of releasable 
vesicles (white patches in Figure IF) whereas the same number 
of 30Hz-stimuli left more release sites functional with one or 
more releasable vesicles (Figure IE). Second, 600 stimuli at 60 Hz 
caused a decrease of the overall number of vesicles per bouton 
to about 84% (mostly red areas in bouton lumen of Figure IF), 
whereas the same number of stimuli applied at 30 Hz resulted 
in a stronger overall decrease to about 73% (mostly yellow in 
Figure IE). Since we modeled vesicle dynamics during synap- 
tic stimulation as an undirected diffusion process that matched 
our experimental observations remarkably well, it seemed likely 
that vesicle diffusion is an important limiting factor during high 
frequency stimulation episodes in vivo (see Discussion). This 
included the disproportionally faster depletion of vesicle release 
at 60 or 80 Hz compared to lower frequencies (Figure 4A). We 
further compared experimental and simulation data in stop and 
go experiments (from Figures 3F,G), where we applied 100 stim- 
uli at the given frequency and 5 s stimulation pause alternatingly 
(see also Figures 3F,G). Data shown in Figures 4B-E showed that 
the peak amplitudes decreased with increasing stimulus number, 
as expected for all stimulation frequencies. We found that both 
experiment and simulation showed a rapid recovery from deple- 
tion during the pausing periods (Figures 4B-E). In both cases the 
recovery was considerably larger during high-frequency stimu- 
lation (Figures 4D,E) than at lower frequencies (Figures 4B,C). 
This could be attributed to the fact that the overall failure rate in 
vesicle release (i.e., a release event is triggered, but no vesicle is 
present at the active zone) was larger at high frequencies, which 
in response causes a slower decline in postsynaptic response. 
This effect was a further indicator for vesicle diffusion limiting 
synaptic response. 

In order to quantify the results in Figure 4A, we applied dif- 
ferent measures (least-squares and L2-norm, see Materials and 
Methods for more detail) to the experimental and simulated data 
(Table 2). These similarities demonstrate that our model includes 
the necessary biophysical laws that explained the stimulation- 
induced vesicular dynamics of multi-release-site boutons in space 
and time and hence might be used to analyze vesicle move- 
ment and their potential restriction for release at high temporal 
and spatial resolution. We saw that the overall bouton dynam- 
ics in vivo and in silico show significant similarities, though 
the diffusion-based in silico results demonstrated a slower decay 
phase at high frequency stimulation in the initial stimulation 
phase, thus indicating that pure diffusion might not suffice for 
high amplitude/high frequency output (see Discussion). 

2.2.2. Bouton performance is differentially affected by bouton size, 
release probability, and stimulation frequency 

Since our three-dimensional model described the vesicle dynam- 
ics in a multi-release-site bouton during example stimulation 



Table 2 | Comparing the experimental and simulated data at different 
frequencies and with the least-squares method S and /.2-norm (see 
Materials and Methods) shows good agreement between 
experimental and simulation data shown in Figure 4A. 



Error norm/frequency 


15 Hz 


30 Hz 


60 Hz 


80 Hz 


Least-squares S 


0.03 


0.02 


0.07 


0.08 


Z_2-norm 


0.11 


0.04 


0.15 


0.1 



frequencies similar to the experimental results, we systemati- 
cally varied each individual component using as standard settings 
P 0 = 6%, d = 3 |im (10 active zones) and co = 40 Hz, with ini- 
tial vesicle density set to 400 (xm -3 , in order to dissect their 
respective roles in determining the functional properties of dif- 
ferent bouton types [Figure 5 (dotted)]. From these simulations 
we derived empirical laws [Figure 5 (lines)] as a heuristic fit of 
the in silico data (see Materials and Methods) in order to assist the 
interpretation of our simulations. 

We found that high P„ -values strongly affected the magni- 
tude of the overall vesicle release at the beginning of a 40 Hz 
stimulation train (Figure 5A). However, the strong initial vesi- 
cle release disappeared within a second of stimulation due to 
the rapid depletion of vesicles at individual release sites. Thus, 
high P 0 's seemed only meaningful if large amplitude outputs 
were needed for short time periods, whereas small P 0 's provided 
long-term reliability of smaller amplitude outputs (Figure 5A). 
The size of boutons affected its performance under otherwise 
constant conditions (P 0 = 6%; co = 40 Hz) on a rather long 
time-scale. Bigger boutons with correspondingly more release 
sites showed larger outputs that depleted significantly slower 
than those of smaller boutons (Figure 5B). Thus, larger bou- 
tons seemed to be better suited for sustained periods of synaptic 
activity whereas smaller boutons with fewer release sites were 
only reliable during short trains. Note that the term "reliable 
transmission" refers to the phases in which vesicle exocytosis 
is constant, thus the phase in which the system is not limited 
by high stimulation frequencies, bouton size or vesicle output 
probability. 

Our analysis of the effects of a broad range of stimulation 
frequencies revealed that frequencies of up to 40-50 Hz were 
reliably transmitted by our standard type lb bouton for sec- 
onds (initial plateaus in Figure 5C). At higher stimulation fre- 
quencies reliable vesicle release was limited to only few stimuli 
until it becomes diffusion-limited (see Figure 5C). Thus, type 
lb boutons could transmit a broad range of stimulus frequen- 
cies even in the absence of vesicle recycling, however, frequen- 
cies above ~50 Hz could only be maintained for very short 
time-periods (Figure 5C). 

The data in Figure 5B, with initial vesicle density set to 
400 |im~ 3 and 7 active zones for the Is bouton, show that the 
exocytosis dynamics are separated into a stable phase in which 
exocytosis is constant and an unstable phase in which there is a 
decay in vesicle exocytosis. The law in Figure 5C shows that in the 
range above 50 Hz the system becomes unstable nearly instanta- 
neously. Experimentally this effect is observed at around 45 Hz for 
an output probability estimated at roughly 7%. 
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FIGURE 5 | Bouton output and output endurance are differentially 
affected by stimulation frequency, bouton size, and P a . Vesicle density is 
set to 400 Lirrr 3 . In each figure row only one of these three parameters were 
varied, while keeping the other two fixed. Systematic analysis of the effects of 
changes in P a (A), bouton diameter (d = 2, 3 and 5 Lim) with 7 10, and 14 
release sites respectively (B) and stimulation frequency (C) on the simulated 
total output dynamics over several simulations (dotted) or on derived 
biophysical depletion laws (lines). If not indicated otherwise we used a P 0 of 



6%, a bouton diameter of 3 Lim and a stimulation frequency of 40 Hz. Note 
that the total bouton output is initially very large at high P D -values, but it rapidly 
depletes to the output of low-P 0 boutons that continue reliable transmission. 
Further note that larger boutons harbor more release sites resulting in larger 
and longer-lasting total outputs (B). Stimulation frequencies up to 40 Hz can 
be transmitted reliably over some time (plateaus in C) while higher 
frequencies can only be maintained for short time periods. Interestingly, the 
long term behavior is independent on the choice of P 0 for P 0 > 5%. 



2.2.3. Decomposition of the compound response of Is/lb-harboring 
NMJs 

Since all three analyzed parameters affected bouton performance 
in distinct manners it was tempting to combine these individual 
model-derived characteristics into realistic functional profiles of 
the two larval bouton types. If one assumed that the active zones 
of the smaller type Is boutons operated with a rather high vesicle 
release probability (in the range of P 0 = 40%) while the bigger 
type lb boutons used a lower, in the range of P 0 = 6%, (Kurdyak 
et al, 1994; Lnenicka and Keshishian, 2000) our above simula- 
tions allowed clear predictions as to how the output of a NMJ 
with roughly similar numbers of both bouton types will look like 
(Figure 6A): a continuous stimulation of, e.g., 40 Hz with initial 
vesicle density set to 400 u.m~ 3 will result in a compound junc- 
tional output (dark gray line) that is initially dominated by the 
strong release from Is boutons (black line) before their contri- 
bution decays due to rapid vesicle depletion at their active zones 
(white patches in Figure 6B). In contrast, the type lb boutons 
of this NMJ are expected to show reliable but small amplitude 
performance (gray line in Figure 6A) with no extensive signs of 
diffusion-based vesicle depletion (left bouton in Figure 6B). 

Intriguingly, in our in vivo experiments in which we stimulated 
the axons of both motor neurons synchronously at various fre- 
quencies (Figure 3E) we observed a very similar behavior of the 
compound response as indicated by an initial steep decrease in 



the eEJP amplitudes followed by an abrupt transition (arrows in 
Figures 3E-G) into a slower decay phase. These results suggested 
both in silico and in vivo that type Is boutons were not well suited 
to reliably transmit sustained high-frequency stimuli and they 
indicate that it was unlikely that in a behaving animal the natural 
motor pattern fires type-Is and -lb boutons synchronously. 

2.2.4. The natural firing pattern matches the functional properties of 
type Is and lb boutons 

Our above results predicted that lb boutons are likely utilized to 
transmit longer-lasting high-frequency stimuli whereas Is bou- 
tons are better suited for few low-frequency events. We tested 
these predictions by analyzing the natural firing pattern that is 
generated by the larval central pattern generator to trigger coor- 
dinated muscle contraction. In order to differentiate between the 
patterns transmitted by Is and lb boutons we established simul- 
taneous recordings from the neighboring muscles 6 and 13 that 
shared a common Is-innervation from the same motor neuron 
but received an independent Ib-innervation each from two dis- 
tinct motor neurons (Lnenicka and Keshishian, 2000). We found 
that the natural motor pattern consists of two discrete compo- 
nents (Figure 6C): a tonic eEJP pattern originating from type-lb 
boutons (arrowheads) that displayed in each train a bell-shaped 
frequency ramp of pulses ranging between 10 and 45 Hz on 
muscle 6 (black line in Figure 6D; 50.8 ± 5.4 pulses per train, 
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FIGURE 6 | Rational design of a naturally performing NMJ in silico. 
(A,B) Simulation of the total (A) or spatially resolved (B) vesicle dynamics 
of a type lb (3|im, 10 active zones; P 0 = 5%) and Is (2|im, 7 active 
zones; P 0 = 40%) bouton during continuous 40 Hz-stimulation. In both 
cases the initial vesicle density was set to 400 |j,m~ 3 . (A) Under these 
conditions the total output dynamics of combined Is and lb boutons 
showed a characteristic biphasic depletion. (B) Is boutons showed more 
depleted release sites than lb boutons at similar stimulus numbers. (C) 
Simultaneous recordings of the natural larval motor pattern in muscles 6 
and 13 identified a synchronous firing pattern (arrows) originating from one 



motor neuron forming Is boutons on both muscles. The remaining 
asynchronous firing (arrowheads) originates from muscle specific motor 
neurons forming lb boutons. (D) Dissected frequency profile of natural 
larval motor patterns plotted as a function of the relative train length. Data 
represent means ± SEM from 10 identified Is- and 10 lb-motor patterns 
from muscle 6. (E) Simulation of the vesicle output of an in silico NMJ (17 
3|im boutons, P 0 = 7%; 19 2|j.m boutons (7 active zones), P 0 = 29%). 
Initial vesicle density was set to 275 liirr 3 . This NMJ exocytosis pattern 
was driven by the brain governed motor pattern as measured 
experimentally and presented in (D). 



mean ± s.e.m. of n = 10). The second component originated 
from type-Is boutons and consists of only few large pulses at 
a rather stable frequency of about 12 Hz (arrows in Figure 6C; 
gray line in Figure 6D; 5 ± 0, 7 pulses per train, n = 10). These 
results demonstrated that Is or lb motor neurons fire with pat- 
terns that corresponded well to our above predictions of their 
individual physiological properties (see Discussion and Figure 7). 
Conversely, when we fed the naturally observed firing pattern of 
both boutons into our model we found that the in silico bou- 
tons reproduced the respective patterns consistently even in the 
absence of vesicle recycling (Figure 6E). 

3. MATERIALS AND METHODS 

3.1. ELECTROPHYSIOLOGICAL RECORDINGS 

Electrophysiological recordings of membrane potentials of body 
wall muscle 6 were performed on filet preparations of size- 
matched mid third instar male larvae at 22C. Dissected larvae 



were first immersed in normal saline (NS) solution containing 
2mM Ca 2 + (] an and Jan, 1976) that was then replaced by NS 
containing 2 \iM Bafilomycin Al (BafAl). After 5 min incubation 
muscle 6 of segment A2-3 was impaled with a 1530 MQ, micro- 
electrode filled with 3 M KC1. For stimulation the cut end of the 
intersegmental nerve of this segment was placed into a suction 
electrode and suprathreshold stimuli were applied at the indi- 
cated frequencies with or without intermitting pausing intervals. 
Cells with a resting potential more negative than —60 mV and an 
input resistance R,„ of 5 MQ, were used for further analysis. Data 
analysis was performed off-line (ClampfitlO, Axon Instruments). 
Recordings of natural motor patterns were performed similarly 
(NS containing 2 mM Ca 2+ without BafAl) except that the motor 
nerves were left attached to the ventral nerve cord and that we 
recorded simultaneously from muscle 6 and muscle 12 or 13. This 
configuration allowed us to differentiate the firing patterns of type 
lb- and 1 s-innervations (Chouhan et al., 2010). Only those trains 
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FIGURE 7 | Parameters defining the configuration of a NMJ. (A) The 

prevailing firing pattern of a neuron sets the standard for the functional 
properties of its synaptic boutons. (B) In order to guarantee reliable 
transmission of a given firing pattern the already existing synaptic boutons 
may modify their genetically determined layout by adapting the here 
identified critical parameters such as the release probability and the bouton 
size (together with the number of release sites). (C) Throughout larval 
development individual boutons are added to the terminals of both motor 
neurons resulting in strings of boutons with typical morphologies. 



were taken for further analysis that correlated with a peristaltic 
contraction wave of the segmentally repeated body wall muscles. 
eEJP-amplitudes in Supplemental Figure SI were measured from 
the extrapolated resting potential to the peak of the eEJP. 

3.2. PARAMETERS USED IN THE MODEL 
3.2.1. Bouton size 

Paraformaldehyde fixed larval filet preparations were immuno- 
labeled with antibodies recognizing a neuronal surface epitope 
(anti-HRP) and the active zone protein Bruchpilot (Wagh et al., 
2006). Confocal images of 12 of such labeled M12/13-NMJs (e.g., 
Figure 1A) were used to estimate the diameters of synaptic bou- 
tons of type Is motor neurons and type lb motor neurons using 
the measurement tool of ImageJ. As exemplified in Figure 1A 
bouton sizes varied more in type lb branches than in type Is 



branches but they were on average larger than type Is bouton. 
Terminal boutons showed particularly large diameters. Because 
of these considerable variances we chose to confine our current 
analysis to three representative classes of bouton sizes with the 
following diameters: Is: 2 \im; lb: 3 |xm; terminal bouton: 5 u,m. 

3.2.2. Vesicle density 

The vesicle densities used in our model were determined as 
follows: (1) The total number of quanta that can be released 
from a Drosophila NMJ has been previously determined to be 
in the range of 85,000 (Delgado et al, 2000). This translates to 
850-1700 vesicles per bouton for NMJs with 50-100 boutons 
and a concentration of 130-270 vesicles/ u,m 3 for an average bou- 
ton (diameter: 2.5 \im, inner vesicle free diameter: 1.5 |xm). (2) 
Ultrathin-sectioned larval boutons showed on average 1-2 vesic- 
ular profiles every 200 nm (Sigrist et al., 2000; Steinert et al., 2006) 
mounting up to about 500-1000 vesicles/|xm 3 . For our simula- 
tion experiments we have chosen a vesicle concentration between 
these estimates: 150-500 vesicles/|im 3 . 

3.2.3. Output probability P 0 

The output probability P 0 values used in this study correspond 
to the experimentally determined vesicle release probabilities P of 
synapses of the two bouton classes Is and lb that were estimated to 
average P values of 29 and 7%, respectively (Kurdyak et al., 1994; 
Cooper et al., 1995). Although a large body of evidence suggests 
that vesicle release probabilities are non-uniform and subject to 
change during facilitation and depression phenomena we chose 
in this stage of our model to simplify and keep the P 0 values con- 
stant during each simulation. To still account for such changes 
we assayed the effects of different P 0 values in a bandwidth of 
1-90%. The P 0 in our model describes the overall likelihood 
that a given vesicle is released from a given active zone and 
therefore represents a summed value for all processes involved 
in an individual vesicle release event including vesicle docking 
and priming, the functional status of voltage gated calcium chan- 
nels, presynaptic calcium dynamics and calcium triggered vesicle 
fusion. 

3.3. DIFFUSION-REACTION MODEL FOR VESICLE DYNAMICS 

The mathematical model presented in this manuscript is based 
on the coarse-scale model of brownian particle motion proposed 
in Einstein (1905). Using this approach we can deduce that vesicle 
dynamics can be described by means of a diffusion-reaction equa- 
tion. To achieve this result, a fixed control region V cR 3 and an 
interval of time X = [0, T] are considered. The total number of 
vesicles contained in V at time t e X is calculated as 

N v (t) = / c(x, t)dv, (1) 
Jv 

where c (units [c] = vesicles • |xm~ 3 ) is the concentration of the 
vesicles. The variation of Ny over time is regulated in principle by 
the flux of vesicles through the boundary 3 V, a source and a sink 
term, i.e., 

d t N v = - I J ndcr+ I (R-A)dv. (2) 
Jav Jv 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



September 2014 | Volume 8 | Article 101 | 9 



Knodel et al. 



Boutons tuned by firing pattern 



Assuming regularity, we can extend / to the interior of the con- 
trol volume V. Using Gauss' Theorem, and a standard localization 
argument, yields the local form of (2): 



d t c : 



-V -J + R-A. 



(3) 



In (2) and (3), / represents the vesicle flux vector, while R and 
A represent, respectively, recycling and absorption of vesicles. In 
general, both _R and A depend on concentration and may also 
exhibit explicit dependence on time and spatial coordinates. Since 
A should vanish when concentration is zero, it is given through 
an expression of the type A(c, x, t) > 0, such that the condi- 
tion A(0, x,t)=0 applies. The physical units of A are [A] = 
vesicles/ (|xm 3 s). For the purposes of our manuscript, recycling 
is switched off (i.e., R is set equal to zero from the outset), and 
only the depletion of vesicles is accounted for. Using Fick's Law, 
and assuming isotropic diffusion, the vesicle flux vector can be 
expressed by / = —DVc, with D the diffusion coefficient, and 
Equation (3) can be rewritten as 



zero Neumann boundary condition for the vesicle flux, the total 
number of vesicles contained in a bouton is not conserved 
over time. To account for the fusion of the vesicles with the 
membrane of the bouton, we introduce a sink term f syn = A 
that, under the conditions discussed below, removes one or no 
vesicle per action potential from each synaptic region of the 
bouton. 

Vesicle exocytosis is a discrete event that reduces the pool of 
vesicles by an integer number of exocytosed vesicles. The sink 
term defined in the following models exocytosis based on the 
general synaptic vesicle release probability P 0 and a stochastic 
probability P, that is calculated for each synapse in each simula- 
tion time step. Since vesicle motion is modeled on the continuous 
scale, and vesicle exocytosis relates to single vesicles, the sink term 
needs to couple the continuous scale with the single particle scale. 
This leads us to the following definition of the vesicle releasing 
sink term/ s} ,„: 



d t c = V • (DVc) - A, Vxe VandVf e T. 



(4) 



fsyn(x, t) :— «o 



Ha i 
J2S(t-t n )h(x,t-)\ 

n=l J 



(8) 



The computational domain, over which (4) is solved, is denoted 
by Q C R 3 and coincides with the region of space occupied by a 
bouton. The set £2 is partitioned in subdomains f2,-, i = 0, . . . , S, 
with S being the number of synapses. £2n is the subdomain 
containing vesicles which are not part of synaptic regions and 
accounts for the largest part of the bouton. For i = 1, 2, ... S, 
the £2; are the synaptic regions. The computational domain Q 
therefore is written as 



i=i 



(5) 



see Figure 2. For all outer boundary surfaces, we use constant 
Neumann zero boundary conditions: 

/ .«= -(DVc) -n = 0, Vx e andVt e 1. (6) 

Finally, we close the mathematical problem by prescribing the 
initial condition 



c(x, 0) = c 0 , Vxe £2, 



where t = 0 denotes the origin of time. 



(7) 



3.4. THE SINK TERM 

The reason we use a sink term rather than a Neumann flux 
boundary condition to simulate vesicle exocytosis is that such 
a boundary condition would describe the situation in which 
the vesicles are permitted to cross the membrane of the bou- 
ton. However, this is not the observed mechanism of vesi- 
cle depletion in a bouton. Indeed, when the signal of the 
action potential reaches a synaptic region, the vesicles avail- 
able for exocytosis fuse with the inner surface of the bouton 
and release neurotransmitters, which, in turn, are allowed to 
cross the membrane of the bouton. Therefore, in spite of a 



h(x, <■-) := ® (Po ~ Pi(t„)) ®+{N Qi (t-) - 1) (9) 
;= l 

Xn,(s)c(s, t~) 

Nq, (tn) 

with ie!2, t n =^el, and n = 1, . . . , Ma- In the 
following, we explain the terms used in the above 
equations: 

1. Given the interval X, that represents the lapse of time during 
which the system is observed, the natural number Ma € N is 
the total number of times (modeled as points of T) in which 
the action potential stimulates the depletion of vesicles. The 
generic instant of time t n e 2, with n = 1, . . . , Ma, corre- 
sponds to the nth stimulus, and is defined by t„ := ^, with 
w being the stimulation frequency. 

2. i refers to the synapse number. 

3. The mathematical definitions of the Dirac Delta S, Heaviside 
distributions 0 and © + , and characteristic function arc 
given in Section 3.9. The Dirac Delta 8 (f — f„) restricts exocy- 
tosis to the time points t n of action potentials; the character- 
istic function xq, restricts the activity of the sink rate solely 
to the ith synaptic region. In the present setting, the Dirac 
Delta and the Heaviside distributions are dimensionless, while 
h(x, t~) has units |xm~ 3 . 

4. P 0 defines the release probability of vesicles at single synapses. 

5. N^it^) — 1 is introduced to ensure that the Heaviside func- 
tion & + is equal to one, if there is at least one vesicle present 
at the corresponding release site. Therefore, in order to obtain 
the result 0+ (Afo,.(f~) - l) = 1, the quantity N Qi (f~) - 1 
has to be non-negative. 

6. A random generator produces values P, between 0 and 1 for 
each synapse at every action potential. If P; is lower than P 0 , 
a vesicle is exocytosed. This is regulated by the expression 
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© (P 0 — Pi(t n )) summed over all synapses i. Thus, in between 
action potentials the model describes a standard diffusion 
process. 

7. Given a generic function of time, (p : I \-+ R, the notation 
<p (t~) stands for <p (t~) = lim f ^ t - <p(t). 

Moreover, the constant an shall be set equal to unity from here 
on. The sink term (8) is now treated in the following way. We 
conventionally set f o = 0 and tj\f A + 1 = T, and write the interval 
J = [0, T] as X = U^L + ! 1 1 m , where 

Ti = [0, h], X q = (t q -i,t q ], q = 2,...,Af A , 

X MA + 1 = (t^ A ,T], (10) 

In each open interval {t m -\, t m ), m = 1, . . . , Na + l,/ syn (x, t) 
vanishes identically and, thus, the solution to Equation (4) rep- 
resents pure diffusion. However, if at a given fp = ^ € I, with 
p € [1, . . . ,Na), the Heaviside functions featured in (9) are 
simultaneously different from zero for the synaptic region £2;, 
the sink term introduces a jump in the value of concentration at 
tp, for all x G £2;. Hence, by employing the notation c(x, tp) = 
lim ± c(x, t), we obtain 



in the interval (f„, f„+ i). We remark that, when the vesicle con- 
centration inside the synaptic regions is such that it is safely 
approximated by its mean value, i.e., N^./VoK^,), for all i = 
I, . . . , S, Equation (11) becomes 



c(x,t+) = c(x,t-) 



1 



Vol(£2,) 

Consequently, Equation (14) can be approximated as follows 

s 

c (x, t+) = c(x,t-)-J2®( p o- Piitn)) 



(15) 



© +(Nn i (t-) - l) 



Xn,-(*) 
Vol(£2,) 



(16) 



We compute the variation of the overall number of vesicles 
contained in the bouton by integrating Equation (14) over £2 
and using the properties of the characteristic functions (see 
Section 3.9), i.e., 



Nfi (*+) = N Q (f,7) - J2 ®( p o ~ Pi(t„))& + {N Ql (t-) - 1). (17) 



(x,t+)=c(x,t-) 



c (x, t p ) 



(11) 



The "updated" concentration c (x, t£) is treated as initial condi- 
tion for the diffusion equation in the interval (tp, fp+ i). 

We remark that, according to (11), only one vesicle is removed 
from the region £2, at tp. Indeed, by integrating (11) over £2,, and 
recalling that (f*) = f Q± c (x, tp) dv, we obtain 



Na, (f+) = No, (t~) 



1. 



(12) 



We denote now by c„, with n = 1, . . . , A^, the restriction of the 
solution to Equation (4) to the interval (t n _i, t n ). By continuat- 
ing each c n to the time point t„, we set c n (x, t n ) = c(x, t~), and 
define the number of vesicles contained in £2,- at time t n in the 
following way: 



/ c n (x, t n )dv = Nn t (t n ). 



(13) 



To evaluate the Heaviside function 0 + of Equation (9), we use 
the number of vesicles NQ^tn), determined in Equation (13). 

More generally, however, for a given x e £2, the jump in the 
vesicle concentration at time f„, with n = 1, . . . , A/a> is given by 

s 

c (x, t+) = c (x, t-) -J> (P 0 - P,(t n )) ® + (N Qi (t-) - 1) 



;= i 



XQ : (x)c(x, t n ) 

Na, (tn) ' 



(14) 



with x G £2. The concentration c(x, f+) of Equation (14) is the 
"initial condition" that has to be used for solving Equation (4) 



We remark that the number of vesicles does not vary, 
i.e., Na(t+)=N n (t-), if either 0(P o -P;(f„)) or 
&+(Na,(t-)- l) is zero. 

In order to determine the solution to Equation (4) in each 
interval I m , with m = 1, . . . , Ma + 1, we consider the restric- 
tion c m to T m , and discretize Equation (4) in time by using a 
first-order finite difference implicit Euler scheme. The result of 
this calculation reads 



c k+l -c k 



fc V • (dv<4 +1 ) , Vx e O, f* t k+l el m (18) 



where c l m (x) = c m (x, t l ), with I = {k, k + 1}, denotes the con- 



centration in x e £2 at time r e X m , and r 



f k+l 



t k > 0 



is the amplitude of the A:th time step. Then, we solve Equation 
(18) by using the Finite Volume Method. For a detailed presen- 
tation of the solution method we refer to standard literature on 
Finite Volumes. In this paper, we sketch very briefly the procedure 
to compute (t~) numerically. To this end, we firstly cover 
the entire computational domain with a finite element triangu- 
lation. The nodes of the triangulation are denoted by x s , with s = 
1 , . . . , 7V„ being the node index, and J\f„ being the total number 
of nodes. The unknown of the problem, i.e., the vesicle concen- 
tration, is approximated by the sum c(x, t) ~ Y^f= \ ^(^^(x), 
where d{t) represents the time-varying value of concentration 
at the node x s , whereas ifr s is referred to as shape function. It 
is such that ir s (x r ) = S s r , with 5* being the Kronecker-delta (it 
returns 1, if r = s, and it returns 0 otherwise). Secondly, we con- 
struct a dual grid by following the methods outlined in Frolkovic 
(1996). An "element" of the dual grid is referred to as finite vol- 
ume. The finite volume associated with the node x s is denoted 
by B 5 in this work, and is obtained by joining the barycenters of 
the elements sharing the node x s with the midpoints of the edges 
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connecting x s with its neighboring nodes. The finite volumes B 5 , 
with s = 1, . . . ,N n , cover the entire computational domain. For 
details the reader is referred to Frolkovic ( 1 996 ) . The whole proce- 
dure is applied also to the synaptic regions £2,-, with i = 1, . . . , S. 
Let then B y - specifically denote the finite volume associated with 
the node x« € £2; of the discretization covering the synaptic 
region £2,, with j = 0, . . . , M; being the nodal index of the dis- 
cretization of £2,, and (M,- + 1) the total number of nodes of 
this discretization. Then, it holds that Vol(£2,) = U.J_ 0 Vol(Bij). 
Moreover, the number of vesicles Nf2 ; (f~) is approximated as 
follows 



Nn,(t n ) = / c n (x,t n )dv 



(19) 



« ^ C n (Xjj, t„ ) V0l(By), Xij £ By C £2;, 

; = 0 

V j = 0, . . . , Mj. 

Finally, by substituting Equation (19) into Equation 
(9), we obtain the discretized form of the expression 

h( Xs ,t-y. 



S , M,- 

h (x,, t~) = Y J & ( p ° - ( X] c„ (x,j, f,;) Vol(B,j) - 1 



xn,-(x s )c n (*s, t„ ) 



(20) 



bouton diameter as d and use characteristic time points with 
i = I , . . . , 6 for tuning the empirical laws. Note that the time- 
points t2, fci, ?6 refer to the time-points in which the respective 
functions reach the value 0.01. 

The characteristic time points originate from heuristic fits 
using logarithmic representations resp. projections in x and y 
space allowing for heuristic determination of approximating 
curves. The same holds true for the E(t, . . . ) which are used in the 
following. All equations stated in the following are used in dimen- 
sionless form, yet for the sake of clarity we omit the procedure of 
making the equations dimensionless during the introduction of 
the following equations. In the following section we define 6 (x) as 



e(x) :-- 



1 x < 0 
0 x > 0 



(21) 



3.7. 1. Bouton size-variation law 

Considering the bouton size to be variable we can investigate the 
effect of the bouton size on the postsynaptic response that was 
simulated in Figure 5C. From the in silico experiments, we derive 
the following empirical law with S = 7, 10, 14 for d = 2, 3, 5 (xm 
respectively and P 0 := 6%: 

E(t, d) = P 0 -S-(6(t- t x {d)) + e(t\{d) - t)) ■ 



with characteristic time points 



3.5. NUMERICAL METHODS 

The model equations were discretized with finite volumes in 
space and an implicit Euler method in time. The resulting sys- 
tems of linear equations were solved with numerical multigrid 
techniques using Gauss-Seidel pre- and post-smoothers and LU- 
decomposition as base solver. The resulting fine grid problems 
consisted of approximately 20,000-40,000 grid nodes and were 
solved on parallel architectures, the HERMIT Supercomputer at 
the HLRS Stuttgart and on the Cekon and Cesari cluster of the 
G-CSC. 



ti(d) = exp log (11.5) + log 



t 2 (d) 



465 
Tl5 



(, , , , /2900\ 
exp (log (80)+ log f — J 



(23) 



(24) 



3.7.2. Frequency-variation law 

The effect of the stimulation frequency on the postsynap- 
tic response was simulated in Figure 5E. From the in silico 
experiments, we derive the following empirical laws with S := 10 
and P 0 := 6% 



3.6. BUILDING THE COMPUTATIONAL DOMAIN 

Bouton morphologies were constructed with Blender and 
ProMesh and discretized with a tetrahedral meshing algorithm 
using TetGen, that is incorporated in ProMesh (Reiter and 
Wittum, 2014). Examples of the 3D geometries can be found in 
Figure 1 . 

3.7. EMPIRICAL LAWS FOR THE IN SILICO RESULTS 

From the simulations performed in Figures 5A,C,E we derived 
fitted empirical laws that reproduce the qualitative features of the 
measured post-synaptic response. The three parameters that were 
varied are bouton size, stimulation frequency and P 0 . In the fol- 
lowing we denote the number of exocytosed vesicle as E(t, w) 
and E(t, d) respectively, the stimulation frequency as a>, the vesi- 
cle output probability as P 0 , the number of synapses as S, the 



E(t, to) = P 0 ■ S- [6(t-t 3 (a)))+0{t 3 (a)) - t) 

log(10P o S) 



exp 



(f-t 3 («))- 



(f 4 (w) - t 3 (w)) 
with frequency-dependent characteristic time points 

t 3 (co) = exp ( - 0. 027(2. 7o - 270)) 
t 4 (o) = exp (8.8 - 0.775 log M) 



(25) 



(26) 
(27) 



3.7.3. P 0 -Variation law 

The effect of the vesicle output probability P„ on the post-synaptic 
response was simulated in Figure 5A. From the in silico experi- 
ments, we derive the following empirical laws with S := 10 and 
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co := 40 Hz 



E(t, P 0 ) 



o(p,-s- c(0) Po-s + e (€(t) - p 0 s) t{i) p„ > 10% 
e(t-t 5 )p 0 -s+9(,t 5 -t)p 0 s 

{t-t 5 (P 0 )) 



exp log (P 0 .S. 100)- 



(f 5 (Po) - t 6 (P„)) 



with the characteristic time points 
t 5 (P 0 ) = 0.3 • exp 
t 6 (P„) = 500 • exp 

where 



1 



(P„ + 0.115) 
2 



(Po • 100) 



1.2 



P 0 < 10% 
(28) 

(29) 
(30) 



£(f) : = 6»(f - 3) (40 exp (- f ° 9 ) + l) + 0(3 - t)9(t - 18) 
•(1.07 • exp ( - 0.25 • (f - 3)) + 1.2) 
+ 9(18 - t)9(t - 50)0.0102 • exp ( - 0.012 ■ t + 5) 

/ log (100) 
+ 0(50 - f)exp I -t ■ & 6Q + 0.24 

3.8. NORMS FOR QUANTIFYING THE SIMULATION RESULTS 

In order to compare the experimental results with the simulation 
results, we applied two norms to the data: 



1. Least-squares S, defined as S := X!"= i 

2. L2-norm, defined as L2 := 



P'-l) , J2^« -Pn-l) 
h<# 



(Pn-Pl) 



(s, — e,) 

where = , with e,- being the experimental data points 

(Si + e,0 

and Sj the simulated data points, p, denote the stimulus number. 

3.9. SUMMARY OF THE MATHEMATICAL EXPRESSIONS USED IN THE 
MODEL 

The symbols S, © and %q { denote, respectively, the Dirac's Delta 
distribution, the Heaviside function and the characteristic func- 
tion associated with the set Q,,. The symbol ®+ is used to modify 
the definition of 0, as explained below. 

Given a set V <Z R d , with R' ; being the d-dimensional 
Euclidean space, the function xv ■ R d - {0, 1} is said to be the 
characteristic function associated with V. It is defined as follows 



xv(x) = { J; fth 



e V, 
otherwise. 



(31) 



Let V C f2 be a subset of Q, and let / : £2 — > R be an integrable 
function over V. Then, the following identity holds true 

/ X v(x)f(x)dv= f f\ V (x)dv, (32) 
Jn Jv 



where f\y represents the restriction of / to V. In particular, if / is 
everywhere equal to unity, then it follows that 



/ xv(x)dv= / dv: 
Jn Jv 



: Vol(V). 



The Heaviside function, 0, is defined as 



0(?) 



1, if t > 0, 
0, if f < 0. 



(33) 



(34) 



If the function is required to return unity also for £ = 0, then the 
definition (34) is extended as follows 



if £ > 0, 
if I < 0. 



(35) 



The Dirac Delta is a distribution. Given an arbitrary function 0 : 
R — >■ R, the Dirac Delta satisfies the important property 



/ 

J A 



8(s-t)cl)(s)ds = 4>(t), 

) 

if <p is continuous at t e R. 



(36) 



4. DISCUSSION 

In this study we have developed a three-dimensional model of 
presynaptic vesicle dynamics to assess the physiological mean- 
ing of experimentally rather inaccessible bouton parameters like 
the bouton size with a given release site density, vesicle diffu- 
sion and the vesicle release probability. In combination with the 
stimulation frequency and an associated release of vesicles these 
parameters represent simplified but fundamental bouton out- 
put functions that should allow a systematic structure/function 
analysis of two very similar classes of glutamatergic boutons. 

4.1. A THREE-DIMENSIONAL MODEL OF LOCAL VESICLE DYNAMICS 

In addition to state of the art modeling approaches to synaptic 
or neuronal functions we used methods that represent the three- 
dimensional morphology and biological processes derived from 
physical principles. Electrical signaling of neurons is often mod- 
eled as a one-dimensional process using the cable equation and 
compartmental one-dimensional space discretization (Hines and 
Carnevale, 2001). There are few approaches in cellular model- 
ing that make use of 2D or 3D modeling techniques (Smart and 
McCammon, 1998; Ross et al., 2000; Meinrenken et al., 2002; Tai 
et al, 2003; Marpeau et al, 2009; Bielecki et al, 2010). Most com- 
monly, stochastic models are employed to address the high com- 
plexity of three-dimensional simulations. A standard approach to 
in silico studies of presynaptic vesicle dynamics is to model the 
transition of vesicles between defined states with preconfigured 
transition-rate constants (Pan and Zucker, 2009), thus reduc- 
ing the three-dimensional process of vesicle motion to a time- 
dependent process (zero-dimensional in space). Considering the 
underlying question of our study presented here, we needed to 
develop a model that takes into account the spatial dynamics of 
vesicles in 3D. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



September 2014 | Volume 8 | Article 101 | 13 



Knodel et al. 



Boutons tuned by firing pattern 



Modeling neurobiological systems requires incorporating the 
physical laws to which these systems abide. The use of determin- 
istic models derived from continuum mechanics theory allows us 
to describe spatio-temporal vesicles dynamics based on physical 
laws which allows for a much more profound understanding of 
the problem compared to the use of heuristic models which are 
afterwards fitted to data. The parameters used for computations 
were taken either from the literature or EM sections and no 
parameter fitting was used. Combined with a three-dimensional 
representation of the morphology using unstructured grids, in 
this case different types of presynaptic boutons, this approach can 
answer questions regarding structure and function relationships. 
The equations are numerically approximated by a finite volumes 
approach and solved on parallel computers. The added mathe- 
matical techniques allowed us to investigate the interplay between 
structure and function at presynaptic terminals in a highly realis- 
tic manner. In addition to the three-dimensional representation 
of the underlying biophysics and three-dimensional morphol- 
ogy, our approach enables us to continually refine the model by 
including more of the systems properties, such as vesicle recycling, 
calcium dynamics and others depending on the current focus. 

4.2. VESICLE DIFFUSION LIMITS THE AVAILABILITY OF PREEXISTING 
VESICLES AT ACTIVE ZONES 

Our simulations modeled vesicle motility during stimulation as a 
simple undirected diffusion process. These simulations matched 
our experimental results (Figures 1E,F, 3, 4) suggesting that vesi- 
cle diffusion describes the motility of preexisting vesicles during 
periods of high frequency stimulation quite accurately. However, 
there is a debate whether or not active and directional transport 
mechanisms based on, e.g., the Myosin/ Actin system are directly 
contributing to synaptic transmission (e.g., Ryan, 1999; Jordan 
et al, 2005; Tokuoka and Goda, 2006; Seabrooke et al., 2010; 
Kisiel et al., 2011; Seabrooke and Stewart, 2011). In particular in 
the Drosophila NMJ system recent evidence has implicated non- 
muscle myosin II and the unconventional myosin VI in mobi- 
lizing synaptic vesicles during synaptic transmission (Seabrooke 
et al, 2010; Kisiel et al, 2011; Seabrooke and Stewart, 2011). In 
our study we focused on high frequency synaptic transmission 
in the absence of vesicle recycling and hence analyzed exclu- 
sively the motility of preexisting mature vesicles during extended 
periods of naturally occurring firing frequencies of 10-80 Hz. 
Under these conditions the modeled slow vesicle diffusion was 
sufficient to rather accurately mimic our eEJP- depletion exper- 
iments (Figures 3, 4) and they showed no obvious requirement 
for the incorporation of a faster vesicular transport mechanism. 
However, our results also showed that vesicle recycling plays a 
major role during prolonged high frequency stimulation as indi- 
cated by the much slower eEJP-depletion kinetics in the absence 
of BafAl (Figure 3D). It seems therefore possible that fast active 
transport mechanisms are primarily involved in vesicle recycling 
processes that rapidly replenish active zones with vesicles during 
periods of prolonged stimulation. In contrast, preexisting vesi- 
cles that may for example be part of the reserve pool of vesicles 
and that have been assayed in this study seem to be recruited to 
active zones primarily by diffusion. Apparently, the glutamater- 
gic boutons of Drosophila NMJs use both active and passive 



vesicular transport mechanisms. Whether and how they interact 
and what their relative contributions to the rather short natu- 
ral activity patterns (Figure 6C) are has to be analyzed in future 
studies. 

Based on our finding that vesicle diffusion is the primary factor 
limiting the replenishment of active zones with preexisting vesi- 
cles, we have assayed its effects on synaptic performance under 
various conditions. Depending on the activity of a given release 
site, which was determined by its release probability, its stimu- 
lation frequency and the typical number of stimuli per pattern, 
the local vesicle concentration dropped rapidly resulting in indi- 
vidual transmission failures and hence in a fast decline of the 
total bouton output (Figures 5A,B,E,F). The size of a bouton 
affects these parameters indirectly by offering, e.g., proportion- 
ally more active zones in larger boutons and simultaneously 
a larger functional reservoir of vesicles resulting in this exam- 
ple in a larger and more enduring bouton output, respectively 
(Figures 5C,D). Reliable transmission of a given activity pat- 
tern therefore depended strongly on the relative settings of these 
parameters. Based on these simulations we predict that a small 
type Is bouton with its smaller number of rather high P 0 -release 
sites is inappropriate for the transmission of long-lasting high 
stimulation frequencies (Figure 5) and should quickly become 
diffusion-limited (white patches in right bouton of Figure 6B). 
This effect is governed by both bouton size and the correspond- 
ing number of release sites. In contrast, a large type lb bouton 
with its larger number of low P 0 -release sites is likely to be 
better adapted to transmit such a pattern (Figure 5; left bou- 
ton in Figure 6B). Indeed, many experimental high frequency 
stimulation paradigms, in which both motor neurons were fired 
simultaneously by a common stimulation electrode showed a 
depression of early compound eEJP- amplitudes before reaching 
plateau values (Delgado et al., 2000). 

Our predictions now revealed that the observed eEJP- 
depression in HFS-trains is likely due to an initial dominance 
and rapid decay of the Is output while the lb output main- 
tains its reliable transmission (Figure 6A). It thus seems obvious 
that in vivo Is and lb motor neurons can not be fired in a 
synchronous manner. Instead, it is more likely that the type-Is 
motor neuron is either fired with a low stimulus frequency or 
with only few stimuli, or both. Type lb motor neurons are, in 
contrast, capable of reliably transmitting longer lasting high fre- 
quency patterns. We found that the natural firing activities of 
both motor neurons (Figures 6C,D) were strikingly consistent 
with the above predictions indicating that motor neuron firing 
patterns and NMJ-bouton properties are well-adapted features of 
this communication pathway. 

4.3. THE PREVAILING FIRING PATTERN SHAPES BOUTON PROPERTIES 

It has been reported in several central and peripheral systems 
including the Drosophila NMJs that long-lasting changes in 
the neuronal firing activity can change the physiological and 
morphological properties of the corresponding output synapses 
(Mollgaard et al., 1971; Hubel et al, 1977; Lnenicka and Atwood, 
1989; Lnenicka et al., 2003; Sigrist et al., 2003). From these obser- 
vations the concept emerges that the prevailing firing pattern of 
a neuron shapes the properties of its output synapses to ensure 
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reliable transmission of its typical patterns. The simulations per- 
formed in this study revealed several tools a neuron can use to 
tune its output performance to its specific needs (Figure 7). 

Under these conditions the developing bouton has to find an 
economic and still effective/reliable compromise between its size, 
i.e., the number of vesicles and release sites and their release prob- 
abilities (Figure 7B). A sustained high-frequency pattern would 
require a relatively large bouton with low release probability 
release sites. In contrast, if the prevailing activity pattern consists 
of only few APs at a rather low frequency an economic realization 
could be a smaller bouton with high-probability release sites. For 
a developing Drosophila NMJ these are the equivalents of type lb 
and type Is boutons, respectively (Figures 1A, 7C). Early in larval 
development, when the innervated muscles are very small only 
few boutons of each type are required to ensure the transmis- 
sion of the motor pattern and an efficient excitation-contraction 
coupling. The larger the muscle grows the less efficient excitation- 
contraction coupling becomes, requiring the proportional addi- 
tion of more boutons of both types (Schuster et al., 1996). In 
addition to this developmental change of boutons, larval NMJs 
can undergo an experience-dependent bouton adaption (Sigrist 
et al., 2003). The morphology of a mature NMJ (Figure 7C) 
therefore likely represents a cumulative readout of its prevailing 
firing pattern throughout development. 

4.4. CONCLUDING REMARKS 

Our three-dimensional functional model of multi-release site 
boutons allowed us to dissect the physiological meaning of exper- 
imentally inaccessible parameters like the bouton size, vesicle dif- 
fusion and the release probability. Based on the simulations with 
the three-dimensional model, we were able to derive an empir- 
ical and analytical description of the complex system subject 
to variations in bouton size, release probability and stimulation 
frequency. Our model sets the framework for a theoretical recon- 
struction and functional dissection of further realistic parameters 
such as vesicle recycling, the membrane potential and its dynam- 
ics during stimulation, the implementation of voltage-dependent 
calcium channels and intracellular calcium transients etc. Many 
of these parameters are experimentally difficult to access. Their 
implementation into a realistic theoretical model will therefore 
enable insights that have previously been impossible. 
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